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ABSTRACT 

We present the first results from two-dimensional simulations of radiatively-efficient accretion of metal-free 
gas onto intermediate-mass black holes. We fix the shape of the spectral energy distribution of the radiation 
produced near the event horizon and study the structure of the irradiated low- angular-momentum accretion 
flow over three orders of magnitude in radius from the black hole, lO^'*- 10^^ cm for a 100 Mq black hole. 
The luminosity of the central source is made to be proportional to the rate at which gas accretes across the 
inner boundary, which we set just inside the sonic radius. We find that photoionization heating and radiation 
pressure modify the structure of the flow. When the ambient gas density is 10^ cm~^, accretion is intermittent 
and on average reduced to 32% of the Eddington-limited rate, two orders of magnitude below the "Bondi" rate 
evaluated ignoring radiation, in agreement with simplified theoretical models. Even if the vicinity of the black 
hole is supplied with high density gas, accretion is rendered inefficient through heating and radiation pressure. 
Subject headings: accretion, accretion disks — black hole physics — galaxies: active — hydrodynamics — 
methods: numerical — quasars: general 



1. INTRODUCTION 

Known stellar-mass and intermediate-mass black hole can- 
didates accrete through mass transfer from a stellar compan- 
ion without exception. Rapidly-accreting massive black hole 
candidates, such as in Seyfert galaxies and quasars, appear 
to be situated close to dusty, molecular, and highly inhomo- 
geneous gas reservoirs, where cold, dense clumps are in free 
fall or are embedded in a hot, ionized, and tenuous confining 
medium. In neither case does the accretion flow resemble the 
laminar flow of the Bondi-Hoyle-Littleton solution. There is 
a distinct possibility, however, that massive black hole pro- 
genitor "seeds," with masses Mbh ^100 Mq (Madau & Rees| 
|2001 ), did pass through an early phase of quasiradial accre- 
tion in protogalaxies, prior to significant dust enrichment and 
the emergence of a pervasive cold phase in the interstellar 
medium (e.g., Volonteri & Rees"200 5| [Alvarez et al.||2006 



2008 i Johnson & Bromm 2007 ; Pelupe ssy et al.|2007[|Di Mat-] 
teo et al.|2008[|Greif et al.|20 08 ). This was the period imme- 
diately following the onset of atomic cooling, the epoch when 
protogalaxies experienced rapid growth due to cosmic infall 
and yet had potential wells too shallow to retain all of their 
newly synthesized metals. 

Since the accretion rate in quasiradial, Bondi-like accretion 
is proportional to gas density, Bondi-like accretion presents an 
appealing formation process for massive black holes: given 
the presence of sufficiently dense gas in a protogalactic nu- 
cleus, Eddington-limited accretion and mass exponentiation 
on a Salpeter time scale would proceed; stellar-size seed black 
holes would reach masses > lO^M© inferred in hi gh-redshift 
quasars in the cosmic time available to them (e.g.,|Haiman & 



|Loeb|200 It [Volonteri et al.|2003trTSiaka & Haiman f2008'). A 
potential flaw in this reasoning was noticed early, namely, that 
in the event of radiatively-efficient accretion, radiative heating 
and radiat ion pressu re may reduce t he accretion rat e ( Shvarts- 
~an|[T9Tl i Buff & McCray|[T974l IHatchett et al.||1976t |Os- 
ti-iker et al. 1976 ; Cowie et al.'1978Vstellingwe rf|1982t|Bep^ 
man 1985 ; Ricotti et al. 2008 ; Alvarez et al. 2008J) The possi- 
bility of a transition to radiative-feedback-driven temporal in- 
termittency, which has been evident in one-dimensional simu- 
lations of Compton-heated accretion onto massive black holes 



in galaxy clusters (| Ciotti & Ostriker|200l| |2007t [Sazonov et 
ar]2005), has further complicated the assessment of quasir- 
adial accretion's viability as a massive black hole formation 
process. If the accretion is indeed intermittent and episodic, 
how long is the duty cycle, how large is the average accre- 
tion rate, and does it permit evolution of seed black holes 
into quasars? Does episodic accretion proceed in bursts, 
where episodes of Eddington-limited rapid accretion alternate 
with near-complete radiative quenching and outflow, as one- 
dimensional models appear to suggest, or does the multidi- 
mensional nature of realistic irradiated accretion flows imply 
a damping of the fluctuations and convergence to a quasi- 
steady state characterized by a reduced accretion rate? 

Here, we present the first results from our program to in- 
vestigate radiatively-efficient accretion onto black holes em- 
bedded in a realistic galactic gaseous medium with the help 
of high-resolution hydrodynamic simulations. We do not at- 
tempt to simulate fluid flow and radiation emission mecha- 
nisms near the event horizon of the black hole; instead, we 
simply fix the shape of the spectral energy distribution of the 
radiation produced near the event horizon and simulate the ir- 
radiated fluid flow on radial scales spanning the sonic radius, 
where the photoheated fluid falls supersonically toward the 
black hole, and a much larger radius in the gas cloud that is 
entirely unaffected by the black hole's radiative output, where 
boundary conditions for the accretion flow are set. This work 
is organized as follows. In §|2]we concisely describe our nu- 
merical algorithm. A more detailed description and a suite of 
tests will be provided in a longer follow-up paper (Couch et 
al., in preparation). In § [3] we describe our results with par- 
ticular focus on the intermittency of the accr etion flow. In § [4 
we compare our results to analytical models ( [Milosavljevic et 
[aL]|2008 , and references therein) and discuss the outlook for 
massive black hole formation through quasiradial accretion. 

2. NUMERICAL ALGORITHM 

Here we describe numerical methodologies that we have 
implemented in the parallel adaptive-mesh-refinement (AMR) 
code FLASH ( [Fryxellet al.|200 0), version 2.5, to facilitate our 
simulations of low-angular momentum accretion flow with a 
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central ionizing source. The hydrodynamic solver used was 
the piecewise parabolic method (PPM), which is a higher- 
order Godunov method. Contact steepening that sharpens 
contact discontinuities in the PPM module was disabled af- 
ter we noticed that it led to spurious flow near stationary ion- 
ization fronts. By default, FLASH does not conserve angular 
momentum in 2D cylindrical simulations. To simulate fluid 
motion in the azimuthal direction, we introduce the specific 
angular momentum / = rv^ as an additional conserved "mass 
scalar" quantity that is advected with the fluid. Equivalence 
with 3D cylindrically symmetric equations of fluid dynamics 
is achieved by adding the centrifugal force to the total accel- 
eration, R = ^grav ~^ ^rad 

2.1. The Central Hole and the Central Ionizing Source 

The simulations were carried out in cylindrical symmetry 
on a rectangular {R^z) grid in two spatial dimensions, where 
R is the distance from the axis of symmetry. In what follows, 
r will denote the distance from the origin, = R^ +z^. A small 
spherical region r < rhoie was excluded from the hydrodynam- 
ical update; we refer to this region as the "hole." A unidirec- 
tional "outflow" boundary condition, allowing only outflow 
from the simulated spatial domain and inflow into the central 
hole, was imposed at r = rhoie- The gravitational acceleration 
was set to that of a Newtonian point mass located at the ori- 
gin of the grid agrav — -GMhoier/r^. The mass of the material 
flowing into the hole was added to the point mass Mhoie- Self 
gravity of the fluid was ignored because, typically, the mass 
of the photoionized gas around the black hole was less than 
1% of the mass of the black hole. 

An isotropic source of ionizing radiation was placed at the 
origin of the grid. The total luminosity of the ionizing source 
between energies /zi^min = 13.6 eV and /zi^max = 100 keV was 
set to be proportional to the accretion rate into the hole. 
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where e is a radiative efficiency which was uniformly set 

to 0.1. The radiative spectral energy distribution was set 

to the power law Lj^ ex and was calculated in 40 

logarithmically- spaced energy bins. 

2.2. Composition, Chemistry, and Thermodynamics 

We simulated a six- species primordial fluid containing hy- 
drogen and helium and their ions: H, H"^, He, He"^, He"^"^, 
and e~. Abundances of the six species were integrated in 
time with a chemical network that includes collisional ion- 
ization of H, He, and He"^, "case B" recombination of H"^, 
He"^, and He"^"^, bidirectional charge exchange between H (H"^) 
and He"^ (He), as well as primary and secondary photoion- 
ization of H and He. Secondary photoionization and heating 
by fast photoel ectrons were implemented in acc ordance with 
the results of ShuU & van Steenberg| ( |1985| ) and |Dalgamo"et 
[aL] ( p^99| ). Photoionization and photoheating rates were cal- 
culated in a photon-conserving manner (see § |2.3| ) from the 
local, extincted radiative fluxes. 

Photoionization heating and collisional cooling of the 
fluid were operator- split from the hydrodynamic update. 
Cooling processes that were taken into account include 
Bremsstrahlung, collisional ionization cooling of H, He, and 
He"^, recombination cooling of H"^, He"^, and He"^"^, dielec- 
tronic recombination cooling of He"^, collisional excitation 
cooling of H and He"^, and Compton cooling to the cosmic 



microwave background. At the high gas densities in the sim- 
ulations, the Compton cooling was negligible. 

Number densities of all six species and the temperature 
of the gas were advanced simultaneously during the chemi- 
cal and thermodynamic update with the Bulirsch-Stoer-type 
semi-implici t ex trapolation mid-point method of Bader &| 
Deuflhard ( 1983| ). Our implementation of the method closely 
follows that provided in the GNU Scientific Library ( |Galassi| 
et al. 2006 ). This method utilizes the Jacobian matrix of the 
partial derivatives of the right-hand side of our seven coupled 
ordinary differential equations. We evaluate the partial deriva- 
tives analytically. Special care was taken to accurately inte- 
grate the abundance of the neutral and partially ionized frac- 
tion in a strongly photoionized gas. We carried the Richard- 
son extrapolation sequence to four points, which was suffi- 
cient for numerical stability. 

In the operator- split alternating hydrodynamical and 
chemothermodynamic update, the hydrodynamic step was 
globally constrained to not exceed 10% of the minimum elec- 
tron abundance change time scale, A^hydro ^ 0.1 ne/\ne\, in 
combination with the usual Courant limit. The chemoth- 
ermodynamic update is carried out in single or in multiple 
and possibly much shorter Bader-Deuflhard steps. The lat- 
ter are individually limited in every cell by the time scale on 
which any of the six abundances and the temperature evolves, 
A^BD ^ minjO.l f7//|^/|,0.1 T/r, A ^hydroj. where / = H, ^~ 
(cf. e.g., |Whalen & Norman|20Q6l ). 



2.3. Line -of -Sight Radiative Transfer and Radiation Pressure 

Here, we discuss our implementation of multifrequency 
photon-conserving line-of-sight radiative transfer and the re- 
sulting radiation pressure force. We calculate the line-of-sight 
optical depth r^^^ from the central source to the faces of any 
rectangular zone along A^ray = 4 separate rays, k = 1, ...,A/ray; 
the rays are equally spaced in angle as seen from the cen- 
tral source. We also calculate the optical depth to absorption 
within the zone along each of the rays fl^\ The effective flux 



at the zone center is calculated via (see, e.g., Mellema et al. 
,2006, and references therein) 



Lu E^rexp(-r(^^^)[l-exp(-ff)] 



47rr^ 



E-^rays 



(2) 



The total photon absorption rates per chemical species / 
within the zone are then calculated via F/ = ai^iyniVFj^/hu, 
where V is the zone volume. In the limit in which the zone is 
very optically thick, f^^^ 1 , we obtain 



^r.>^(M-iF,(exp(-r^^))^], 



(3) 



where ft is the angular size of the zone as seen from the central 
source, and the extinction external to the zone is averaged over 
rays. The approximation in equation ^ becomes exact in the 



limit Nr 



rays 



oo, as it has to be to conserve photons. 



We calculate radiation pressure associated with photoion- 
ization and Thomson scattering. The radiation pressure force 
is added to the gravitational force of the central point mass. 
We ignore line resonance radiation altogeth er and also ignor e 
any diffusion of ionizing radiation (see, e.g., |Ritzerveld|2005] ) . 



2.4. Mesh Refinement and Initial Conditions 

We simulate a cylindrical domain in the upper hemisphere 
< (/?,z) < 1 pc with 18 levels of mesh refinement. Since 



MILOSAVLJEVIC, COUCH, & BROMM 



3 




Fig. 1. — Central accretion rate and radiative luminosity as a function of 
time for radiative efficiency e = 0.1. 



each AMR block contains 8x8 zones, the maximum spatial 
resolution in the simulation is Ax ^3x10^^ cm. We require 
that the central hole, with radius rhoie = 10^"^ cm, is always re- 
solved at the highest level of mesh refinement. We achieve 
pseudo-logarithmic gridding by capping the resolution at ra- 
dius r with Ax > ^rjr where we choose 7/ = 0.1; this prevents 
use of excessive resolution far from the central hole. The sim- 
ulation domain initially contained uniform-density partially 
ionized gas, with initial electron abundance of = 0-5, at 
temperature 10"^ K and density = 10^ cm~^. The initial 
value for the central point mass was 100 Mq. 

3. RESULTS 

The central accretion rate, shown in Figure [T] oscillates 
between values close to, and occasionally mildly exceed- 
ing the Eddington limit, Mmax ~ ^Edd ~ 2 x 1O~^M0 yr~\ 
and a rate lower by an order of magnitude, Mmin ~ 2 x 
1O~^M0 yr~^ The accretion is approximately periodic with a 
mild trend toward an increasing period separating consecutive 
peak episodes; the period varies between 250 yr and 350 yr. 
Accretion rate falloff following a maximum is roughly expo- 
nential, as is expected since photoionization radiation pres- 
sure in the ionized region surrounding the black hole evac- 
uates the accretin g gas and drives do wn the accretion rate 
(see § 3 in Milosa v^evic et al^|2008[ ). The average accre- 
tion rate and luminosity are (M) = 4.6 x 10^^ g cm~^ and (L) = 
4.2 X 10^^ erg s~\ which says that on average, the black hole 
accretes at 32% of the Eddington limit. The average accretion 
rate is still only ~ 1% of the adiabatic "Bondi" accretion rate 
^Bondi = 7r(GM\y\^)'^nmp/cl(oo), calculated ignoring radiative 
feedback, for an ambient sound speed of Cs(oo) = 14 km s~^ 

During a central accretion maximum, photoionization radi- 
ation pressure drives an outflow in the ionized gas within the 
H II region that has neutral fractions xh ~ 10""^ - 10~^ (Fig-|2l 
lower panels). This leads to rarefaction and exponential drop 
in central accretion. Meanwhile, as radiation pressure sub- 
sides, gas near the edge of the H II region accelerates inward. 
This acceleration is driven by a gas pressure imbalance near 
the edge; the imbalance was inherited from the preceding ac- 
cretion maximum when an outward radiation pressure force 
balanced an inward gas pressure gradient force. The outflow 
intersects with the inflow, and the inflowing gas ultimately ar- 
rives at the edge of the central hole and gives rise to a new 
accretion maximum. The longest timescale in the cycle is the 
inward acceleration time, which is here a factor of ~ 3 times 
shorter than the radial sound crossing time of the H II region 
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Fig. 2. — Gas number density n (upper panels) and neutral fraction xh 
(lower panels) at two separate instances around f ~ 4, 000 yr. The left panels 
show the flow during a central accretion minimum; infalling gas is clearly 
visible at r ~ 1.5 x 10^^ cm. The right panels show the flow during a central 
accretion maximum. The structure near the central axis, in the conical region 
marked by dashed lines, is an artifact of our having set the optical depth to 
zero along the axis, for < tan 5°, to avoid numerical issues arising from 
the coordinate singularity. 



with typical temperature The ~ 4 x 10"^ K. 

Figure ^ shows that the radius of the H II region sur- 
rounding the black hole reaches an asymptotic value of rion ^ 
6 X 10^^ cm. This only somewhat smaller, by 40%, than the 
crude estimate obtained for hydro gen-only cloud in equation 
(4) of jMilosavljevic tCsI] ( [2008| ). The region r < rion is not 
entirely free of neutral gas; during accretion minima, as seen 
in the two left panels of Figure [2j dense clumps form in the 
gas that has been expelled by radiation pressure and is now 
falling back toward the black hole, and the regions in the 
dense clumps' shadows recombine. The range of radii that is 
shaded in Figure|3]2 contains a multiphase medium. The dens- 
est neutral clumps that form near the edges of the H II region 
are six times denser than the ambient gas, but this density is 
not sufficient to resist photoevaporation during accretion max- 
ima. Figure[3|? shows that the flow into the central hole is typi- 
cally supersonic with sonic radius located at Ts ~ 2 x 10^^ cm, 
but during accretion maxima, shock heating in the infalling 
gas drives subsonic flow to the brink of the excised hole at 
^hoie = 10^"^ cm. 

4. CONCLUSIONS AND DISCUSSION 

We have simulated accretion from a uniform, high- 
density (10^ cm~^)^ metal-free protogalactic cloud onto an 
intermediate-mass black hole. At radii from the black hole 
smaller than the size of the excised region, 10^"^ cm, the ac- 
cretion was assumed to be radiatively efficient. We find that 
a compact H II region forms around the black hole. The ac- 
cretion is intermittent due to alternating photoionization radi- 
ation pressure-driven expulsion and external pressure-driven 

^ Atomic densities have been found to reach comparable high levels in 
simulation s of protogalactic collapse (e.g., |Bromm & Loeb|2003| |Wise &| 



Abel 2007l[Wise et al.|2008) . 
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Fig. 3. — Minimum radius at which neutral gas, with ionization fraction 
Xh+ < 5 is found (lower boundary of the shaded region in panel [a]) and 
maximum radius at which ionized gas, with xh+ > 5 ' found (upper bound- 
ary of the shaded region). Panel (b) shows the minimum radius at which the 
radial inflow is subsonic. Because of transient shock heating in the flow, the 



sonic radius, which normally resides at rg ' 
edge of the hole at rhoie = 10^^ cm. 



' 2.2 X 10^"^ cm, shrinks to the 



fallback. The average accretion rate is < 1% of the Bondi ac- 
cretion rate calculated ignoring the radiation's influence, and 
~ I of the Eddington-limited rate. 

The numerical results are consistent with the predictions of 
our toy m odel for episodic quasiradial accretion (Milosavl-J 
[jevic et al.| |20Q8), where we suggested that photoioniza- 
tion heating and photoionization radiation pressure within the 
compact H II region are the dominant accretion-suppression 
mechanisms. We also suggested that Lya line resonance ra- 
diation pressure, which we do not account for in the simula- 
tions presented here, should further reduce the accretion rate. 
The intermittency is qualitatively similar to that seen in one- 
dimensional simulations of the accretion of hot irradiated in- 
tergalactic medium onto black holes in central cluster galaxies 



by [Sazonov eTaTI ( [2005] ) and |Ciotti & Ostriker] ( [2007] ), with 
the notable caveat that photoionization radiation pressure was 
unimportant in that context. 

The simulation in which the accreting fluid was endowed 
with small angular momentum corresponding to maximum 
circularization radius equal to rhoie exhibited similar intermit- 
tency and o nly a slightly reduced a verage accretion rate, con- 
sistent with [Krumholz et al.| ( [2005 , and references therein). 
In contrast with the nonrotating simulations that did not con- 
tain axial outflows, narrow polar outflows form that resemble 
some of those found by Proga ( 2007 ) and Pr oga et al. ( 2008] ) 
in simulations of accretion flows irradiated by a quasar. The 
accretion intermittency and inflow-outflow behavior that we 
observe is qualitatively similar to that seen in Proga et al., but 
additional tests are needed to ascertain whether out axial out- 
flows are physical. 

The overall goal is to combine the highly-resolved simu- 
lations presented here with large-scale cosmological simula- 
tions of how the first galaxie s formed at the end of the cos- 
mic dark ages (e.g., Greif et al. 2008 ). Those simulations can 
resolve the central gas flows to within ~ 0.01 pc, where we 
can implement sub-grid prescriptions based on the calcula- 
tions done here. An important challenge will be to derive pre- 
scriptions for the feedback from black hole accretion that are 
applicable to a wide range of protogalactic conditions. We 
will report on these efforts elsewhere (Couch et al., in prepa- 
ration). 
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